AD-A220  416 


* 


\ 


HARVARD  UNIVERSITY 


Department  of  Statistics 


The  Trick ett- Welch  “Solution*  to  the  Behrens-Fisher 
Problem  Applied  to  One-sided  Tolerance  Limits 
for  Random  Effects  Models 


Mark  G.  Vangel 


U.S.  Army  Materials  Technology  Laboratory 
Watertown,  MA  02171-0001 

and 


DTIC 


F.LKCTE 
APR  12  1990 


Department  of  Statistics 
Harvard  University 
Cambridge,  MA  02138 


Technical  Report  No.  ONR-C-5 
January  1990 


Reproduction  in  whole  or  in  part  is  permitted  for  any 
purpose  of  the  United  States  Government 

This  document  has  been  approved  for  public  release  and 
sale,  its  distribution  is  unlimited. 


90  o v  n  no 


I 


Abstract 

Let  X  be  a  normally  distributed  random  variable  with  mean  y  and  variance 
or2  =  of  4-  a*  .  A  lower  confidence  limit  for  a  quantile  of  this  population 
(i.e.,  a  tolerance  limit)  is  to  be  determined  using  data  from  a  one-way  bal¬ 
anced  random  effects  ANOVA  sample  with  between-group  and  within-group 
variances  of  and  of  respectively. 

For  example,  let  X  represent  the  strength  of  a  randomly  selected  spec¬ 
imen  of  a  material  manufactured  in  a  batch  which  can  be  considered  to 
be  randomly  selected  from  a  population  of  batches.  A  quantity  of  inter¬ 
est  to  aircraft  designers  is  the  ‘B-basis  value’,  which  is  a  95  percent  lower 
confidence  limit  on  the  tenth  percentile  of  the  distribution  of  X.  For  this 
situation,  it  is  important  that  nearly  the  nominal  coverage  probability  be 
attained  whatever  the  unknown  population  variance  ratio.  It  is  also  very 
desirable  that  the  calculated  limit  be  as  large  as  possible,  since  unnecessarily 
low  values  cause  undue  conservatism  in  design. 

This  problem  is  closely  related  to  the  Behrens-Fisher  problem.  We  have 
in  the  one-way  ANOVA  two  independent  mean  squares  with  expected  values 
equal  to  linear  combinations  of  the  variance  components,  while  for  the  two 
sample  problem  the  expected  values  are  the  sample  variances.  The  Welch 
series  approach  (Welch,  1947)  can  be  applied  here  to  produce  an  approxima¬ 
tion  which  is  adequate  for  many  batches.  A  little  known  paper  by  Trickett 
and  Welch  (1954)  describes  an  equivalent  integral  equation  which  is  applied 
to  the  tolerance  limit  problem  with  dramatic  results.  Unlike  the  Welch  series 
calculations,  which  axe  as  tedious  to  do  today  as  they  were  forty  years  ago, 
the  Trickett- Welch  approach  is  numerical,  and  one  can  calculate  the  succes¬ 
sive  approximations  to  orders  inconceivable  in  1954.  In  fact,  though  there  is 
strictly  speaking  no  ‘well  behaved’  exact  solution  to  this  problem,  one  can 
get  amazihgly  close  to  the  nominal  coverage  probability  for  any  value  of  the 
nuisance  parameter  by  beginning  with  the  first  order  Welch  approximation 
and  iterating  an  improvement  of  the  Trickett- Welch  procedure  numerically. 

A  solution  due  to  Mee  and  Owen  (1983)  based  on  Satterthwaite’s  (1946) 
approximation  for  the  distribution  of  a  linear  combination  of  x7  random 
variables  is  compared  with  the  above  approach  and  with  the  solution  for 
known  variance  ratio.  The  comparison  is  made  both  in  terms  of  the  cov¬ 
erage  probability  and  by  calculation  of  the  probability  distributions  of  the 
tolerance  limits. 


1  Introduction 


If  a  material  is  manufactured  in  many  large  batches  and  the  population  of 
interest  consists  of  all  batches,  a  random  effects  model  may  be  an  appropriate 
model  for  measurements  made  on  characteristics  of  the  material. 

Let  Xij  denote  the  jth  of  J  observations  from  the  tth  of  I  batches.  If 
Xij  follows  a  one-way  balanced  random-effects  model,  then 

Xij  =  P  +  b,  +  ,  (1) 

where  p  denotes  the  population  mean,  p  +  6,-  denotes  the  mean  of  the  ith 
batch,  and  e,j  is  the  error  term.  The  6;’ s  and  the  etj’s  are  assumed  to 
be  independently  distributed  normal  with  mean  zero  and  variance  a\  and 
a\  respectively.  An  observation  X  from  this  population  is  thus  normally 
distributed  with  mean  p  and  variance 

ax  —  ab  +  ae-  (2) 

This  paper  presents  techniques  for  determining  one-sided  tolerance  limits 
for  X  based  on  a  random  sample  of  J  items  from  each  of  I  batches.  A  (/?,  7) 
lower  tolerance  limit  is  a  random  variable  T  such  that  at  least  a  proportion  (3 
of  the  population  is  covered  by  the  interval  ( T ,  00)  with  probability  at  least 
7.  The  methods  developed  here  for  lower  tolerance  limits  can  be  adapted  in 
an  obvious  way  to  upper  limits.  We  will  refer  to  /?  as  the  coverage  and  7  as 
the  coverage  probability. 

An  important  industrial  application  of  tolerance  limits  is  to  the  charac¬ 
terization  and  certification  of  structural  materials  for  aircraft.  In  order  to 
determine  the  acceptability  of  material  for  aircraft  applications,  designers 
use  ‘material  basis  properties’  which  are  tolerance  limits  on  the  strength 
of  a  material  as  determined  from  experimental  failure  data.  A  (.90,  .95) 
lower  tolerance  limit  is  called  a  ‘B-basis’  value  or  ‘B-allowable’.  The  more 
stringent  (.99,  .95)  limit  is  referred  to  as  an  ‘A-basis’  value  or  ‘A-allowable’. 

There  is  increasing  interest  in  the  use  of  composite  materials  as 
lightweight  alternatives  to  metals  for  aircraft  applications.  Composite  ma¬ 
terial  properties  typically  exhibit  far  more  batch- to- batch  variability  than 
do  metals;  consequently  there  is  a  growing  need  for  methods  to  determine 
one  sided  tolerance  limits  in  the  presence  of  batch-to-batch  variation. 

Various  approaches  to  this  tolerance  limit  problem  will  be  discussed  be¬ 
low.  The  ratio  of  population  variance  components  is  a  nuisance  parameter 
for  this  problem.  How  one  chooses  to  address  this  complication  is  a  distin¬ 
guishing  feature  of  the  alternative  methods. 
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Section  2  presents  a  method  due  to  Mee  and  Owen  (1983)  based  on 
Satterthwaite’s  (1946)  approximate  distribution  for  a  linear  combination  of 
mean  squares.  A  modification  is  suggested  which  slightly  improves  on  the 
Mee- Owen  result. 

In  Section  3  the  exact  solution  is  derived  for  known  variance  ratio  in 
terms  of  a  generalization  of  the  noncentral  t-distribution.  This  distribution 
is  used  in  Section  4  to  examine  the  effect  on  the  coverage  probability  7  of 
pooling  and  using  a  simple  random  sample  procedure  when  the  variance  ratio 
is  not  zero.  Section  5  consists  of  an  asymptotic  series  solution  (following 
Welch  (1947))  for  the  tolerance  limit  limit  factor  to  terms  of  0(l/n).  In 
Section  6  this  problem  is  formulated  as  an  integral  equation  and  a  method 
due  to  Trickett  and  Welch  (1954)  is  applied.  Following  the  improvements  of 
Section  7,  this  approach  is  shown  to  provide  virtually  exact  tolerance  limits 
for  the  case  of  an  unknown  nuisance  parameter. 

The  Trickett- Welch  approach  has  received  little  attention  in  the  statistics 
literature.  The  dramatic  success  of  this  numerical  method  for  the  problem 
considered  in  this  paper  suggests  that  one  might  profitably  apply  this  tech¬ 
nique  to  a  large  class  of  inference  problems.  Two  such  potential  applications 
are  outlined  in  Section  8. 

The  cumulative  distribution  function  of  the  tolerance  limits  are  examined 
in  Section  9.  These  cdfs  may  be  easily  calculated  in  terms  of  the  generalized 
noncentral  t-distribution  of  Section  3. 

The  discussion  of  the  various  tolerance  limit  procedures  in  Section  10 
makes  use  of  both  the  calculated  coverage  probability  as  a  function  of  the 
nuisance  parameter  and  the  cdfs  of  the  tolerance  limits  in  making  compar¬ 
isons. 

Finally,  a  simulated  dat.'*  r  '  mple  is  considered  in  Section  11. 


2  The  Mee-Owen  Procedure 

Let  n  =  IJ  denote  the  sample  size.  The  parameters  n,  of  and  of  of  the 
random  effects  model  may  be  estimated  by  the  pooled  mean  /t,  the  within 
batch  mean  square  M Se,  and  a  linear  combination  of  MSe  with  the  between 
batch  mean  square  MSb  where: 


I  J  Y 
1=1 ;=1 


(3) 


2 


(4) 


»aj 


and 


MS  ~  ff  ^  ~X-^2 

'  hh  ^-1)  ' 

An  unbiased  estimator  of  the  population  variance  a\  is 
b\  =  MSb/J  +  (1  -  1  /J)MSe. 


(5) 

(6) 


For  0  <  0  <  1,  let  Kp  be  the  (3  quantile  of  the  standard  normal  distri¬ 
bution,  i.e 


(7) 


A  (/3,7)  lower  tolerance  limit  is  a  IOO7  percent  lower  confidence  bound  for 


fi-Kpcrx-  (8) 

By  analogy  with  the  single  sample  case  (see,  for  example,  Owen  (1968)),  one 
seeks  an  estimator  of  the  form 


M  ~  k&x. 


where  k  is  chosen  to  satisfy 

P{p.  -  kbx  <  M  “  K0ox)  =  7- 
Since  fi  is  distributed  normal  with  mean  n  and  variance 

=  (Jo*  +<??)/", 


one  may  rewrite  (  10)  as 


where 


B 


I  R  +  l 

■  V  JR  +  1’ 


(9) 

(10) 

(11) 

(12) 

(13) 

(14) 
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and 

R=.alla2t.  (15) 

The  random  variable  is  approximately  distributed  as  the  ratio 

of  a  x2  to  its  degrees  of  freedom,  where  the  degrees  of  freedom  are  given  by 
Satterthwaite’s  (1946)  approximation: 


/  = 


(*+l)2 


,  1-1 IJ 


l-T 


(16) 


If  7T1( 7,<5)  denotes  the  inverse  of  the  noncentral  t-distribution  with  /  de¬ 
grees  of  freedom  and  noncentrality  parameter  6  then  one  may  make  the 
following  approximation: 

T7*(7,  y/iiKpB) 

k  %  -i - =— - 

y/nB 

Unfortunately,  this  tolerance  limit  factor  k  depends  on  the  nuisance  param 
eter  R.  Mee  and  Owen  (1983)  suggest  replacing  R  with 


(18) 


where  Fn  is  the  100 77  percentile  of  an  F  random  variable  with  numerator 
and  denominator  degrees  of  freedom  /( J  -  1)  and  /  —  1,  respectively.  Rr,  is 
a  1007/  percent  upper  confidence  bound  estimate  for  R  (Searle,  1971,  p.414). 

Having  made  the  approximation  (  17),  we  may  determine  the  coverage 
probability 


P(fi  -  k(R,,)&x  <  H  -  Kpox)  =  lm{r},I,J,R).  (19) 


As  R  tends  to  infinity,  (  17)  becomes 

k  —  Tf2i{’l,'/lKp)/y/I  (20) 

for  all  7/  and  J.  The  case  of  infinite  R  corresponds  to  a\  =  0;  the  model  (  1) 
reduces  to  a  simple  random  sample  of  size  /,  and  (  20)  provides  an  exact 
solution.  Hence  for  all  7/,  /,  and  J 


Jim  J,R)  =  7-  (21) 

n— *00 

For  R  sufficiently  small,  it  can  be  shown  that  7*  exceeds  the  nominal  level 
7.  However,  in  general  7*  will  be  less  than  7  for  an  interval  of  intermediate 
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R  values.  It  turns  out  that  one  can  determine  jj  =  rj(0, 7)  numerically  so 
that  7*  >  7  for  all  /,  J,  and  R.  These  t?  values  are  reproduced  from  Mee 
and  Owen  (1983)  for  various  combinations  of  0  and  7  in  Table  1. 

The  first  improvement  to  this  tolerance  limit  procedure  that  we  will 
consider  is  to  allow  r)  to  vary  with  I  and  J.  This  gives  a  modified  Mee- 
Owen  procedure  with  coverage  probability  closer  to  the  nominal  value.  The 
result  of  this  numerical  work  for  the  case  of  0  =  .90  and  7  =  .95  is  presented 
in  Table  2. 


3  An  Exact  Solution  for  Known  R 


For  a  simple  random  sample,  a  solution  to  the  one  sided  tolerance  limit 
problem  is  readily  obtained  in  terms  of  the  noncentral  t-distribution.  If  one 
assumes  that  the  variance  ratio,  R,  is  known,  then  the  corresponding  prob¬ 
lem  for  a  sample  from  a  balanced  random  effects  model  can  be  solved  almost 
as  easily.  What  is  required  is  the  distribution  of  a  ‘generalized  noncentral  t’ 
random  variable,  a  generalization  of  the  noncentral  t  to  a  random  variable 
with  the  square  root  of  a  linear  combination  of  two  x2’s  in  the  denominator. 

Let  6  =  •JnKpB,  nx  -  I  -  1,  and  ni  =  I{J  -  1)  where  B  is  defined  in 
(  14).  If  R  is  known,  the  tolerance  limit  factor  k  is  the  appropriate  quantile 
of  the  distribution  of 


A  =  (ni  +  n2)1/2 


Z  +  6 

s/  d\Y\  +  d2Y2 


(22) 


where  Z  has  a  standard  normal  distribution;  Yi  is  distributed  as  a  x2  with 
nt  degrees  of  freedom  for  i  =  1,2;  d\,  d2,  and  <5  are  constants  with  di  and  d2 
positive;  and  Z,  Fi,  and  Y2  are  mutually  independent.  Once  this  distribution 
has  been  determined  the  tolerance  limit  may  be  obtained  exactly.  The  cdf 
of  the  linear  combination  Y  =  d\Y\  +  d2Y2  is  show  in  Fleiss  (1971)  to  be 

My)  =  Ev  [X’1+nj  {yf(dxX  +  d2(l~  *)))] ,  (23) 


where  Xni+nj  *s  t^e  chi-square  cumulative  distribution  with  /  degrees  of 
freedom  and  the  expectation  is  with  respect  to  a  beta  random  variable,  X, 
with  parameter  v  —  (n\/2,n2/2). 

By  conditioning  on  the  denominator  of  (  22)  one  sees  that 


FA(k)  =  P(A  <k)  =  E, 


-^"ni+na  (*^,*  +  <*2(1 -*),«)],  (24) 
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where  Tj(t,6)  denotes  the  noncentral  t  cumulative  distribution  with  /  de¬ 
grees  of  freedom  and  noncentrality  parameter  6,  i.e. 

Tj(t,6)  =  jH  $  (t^  -  <5)  Cf  ( w)dw ,  (25) 


where  C/  denotes  the  \2  density  with  /  degrees  of  freedom  and  $(•)  is  the 
standard  normal  distribution. 

For  the  tolerance  limit  problem,  let 


(nj  +n2)I 

d'=  '  rrr 


(26) 


and 


d2 


ni  +  n-i 

jr- 1-  r 


(27) 


where  /,  J,  Kp,  and  R  are  as  in  Sections  1  and  2  and  nj,  ri2,  and  6  are 
given  above.  The  value  k(R)  such  that  Fy\(k)  =  7  thus  provides  an  exact 
solution  to  the  problem  of  known  R. 

Although  the  above  derivation  is  simple,  it  is  apparently  not  well  known. 
A  much  more  complicated  representation  of  the  distribution  of  the  random 
variable  (  22)  is  developed  in  Ray  and  Pitman  (1961). 


4  The  Effect  of  Pooling  on  the  Coverage 
Probability 

The  tolerance  limit  procedure  discussed  in  Section  2  is  conservative  (i.e. 
provides  a  coverage  probability  greater  than  the  nominal  value)  when  the 
population  variance  ratio,  R,  is  small.  Mee  and  Owen  (1983)  therefore  sug¬ 
gest  that  data  be  pooled  and  that  a  single  sample  method  be  applied  when 
the  mean  square  ratio  is  less  than  one.  They  then  proceed  to  investigate  the 
conditional  behavior  of  their  proposed  estimator. 

Using  the  distribution  developed  in  Section  3,  we  shall  determine  the 
coverage  probability  for  a  single  sample  procedure  applied  to  pooled  data  as 
a  function  of  the  variance  ratio.  This  result  will  be  used  to  determine  the 
unconditional  coverage  probability  of  the  Mee-Owen  method  in  Section  10. 

Let  Vj  and  Y2  be  as  in  (  22)  and  let  nj  and  712  denote  the  between  and 
within  batch  degrees  of  freedom  respectively.  The  pooled  variance  estimate 

,2  E/=i  E/=i(*u  -  A)2 

5p  -  —x  ’  <28) 
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where  n  denotes  the  pooled  sample  size,  I  the  number  of  batches,  J  the 
batch  size,  and  p  the  grand  mean.  Partitioning  the  total  mean  square  and 
substituting  (  11)  for  the  variance  of  p,  one  obtains 

Sp _ n  a\Yi  -f  (J&1  +  a\)Yx 

crp}  n  -  1  Jo%  +  cr} 

If  fco  denotes  the  single  sample  tolerance  limit  factor  (e.g.  Owen,  1968,  pp. 
446-448),  then  the  coverage  probability  as  a  function  of  R  is 


7 (R)  =  P(A  -  kQSp  <  p  -  Kaox ) 


with  notation  as  in  Section  2.  Substituting  (  29)  into  (  30)  and  employing 
the  distribution  (  24),  one  may  readily  examine  7 (R)  numerically.  For  the 
present  application  of  (  24)  the  constants  dj  and  can  easily  be  determined 
from  inspection  of  (  30);  they  are  different  than  the  values  assigned  in  (  26) 
and  (  27). 

From  the  typical  plot  in  Figure  1  it  is  apparent  that  the  coverage  prob¬ 
ability  obtained  will  be  substantially  less  than  the  nominal  value  even  for 
fairly  small  values  of  R.  Clearly,  criteria  which  result  in  the  decision  to 
pool  must  be  considered  carefully  if  one  is  to  minimize  the  risk  of  very  an- 
ticoTis'"”"\tive  tolerance  limits  in  *he  presence  of  batch-to-batch  variation. 
Alternatively,  one  might  seek  an  estimator  which  performs  well  for  all  R, 
eliminating  the  consideration  of  pooling  altogether.  This  approach  will  be 
taken  in  Section  6. 


5  The  Solution  for  Unknown  R :  Welch- Aspin 
Series 

For  unknown  variance  ratio,  the  tolerance  limit  problem  is  closely  related  to 
the  Behrens-Fisher  problem.  Since  it  has  been  shown  by  Linnik  (1968,  Ch. 
9  )  that  there  is  no  ‘well  behaved’  solution  to  the  Behrens-Fisher  problem, 
it  follows  that  we  also  are  faced  with  a  problem  without  an  exact  solution. 
However,  one  can  proceed  as  if  a  solution  does  exist  and  attempt  to  approxi¬ 
mate  it.  Following  the  work  of  Welch  (1947)  and  Trickett  and  Welch  (1954), 
two  forms  for  such  an  approximate  solution  are  obtained. 
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A  series  solution  is  developed  first.  While  computationally  simple,  the 
first  order  approximation  presented  here  is  anticonservative  and  may  only 
be  suitable  for  many  batches. 

One  could  improve  this  procedure  by  taking  higher  order  approxima¬ 
tions.  However,  this  becomes  very  tedious  to  carry  out.  Alternatively,  the 
tolerance  limit  factor  as  a  function  of  the  mean  square  ratio  may  be  ob¬ 
tained  approximately  as  the  solution  of  an  integral  equation.  Although  this 
requires  the  use  of  a  computer,  the  method  which  results  appears  to  give 
very  nearly  the  nominal  coverage  probability  -  even  for  small  sample  sizes. 

To  simplify  the  notation  in  what  follows,  let  Sf  be  the  mean  squares,  <r2 
their  expected  values,  and  n,  the  asociated  degrees  of  freedom  for  i  =  1,2, 

•Sf  =  MSfc,  of  =  Jerf  +  erf,  ns  =7-1, 

Sf  =  MSe,  of  =  erf,  n2  =  I(J  -  1). 

The  pooled  sample  size  is  n  —  IJ  and  the  population  variance  is  denoted  by 
cr2  =  a\  =  of  +  of  -  a\jj  +  cr|(l  -  1/J),  (31) 

and  estimated  by 

S2  =  b\  =  Sf/J  +  Sf(l  -  l/J).  (32) 

The  subscript  X  for  the  population  variance  and  estimates  of  this  variance 
will  be  omitted  for  the  remainder  of  this  section. 

The  tolerance  limit  factor  will  be  denoted  k  as  in  (  10)  ,  and  we  define 
A(Sf,S2)  to  be  kb.  The  tolerance  limit  may  be  expressed  as  an  expectation 
with  respect  to  the  distributions  of  the  mean  squares  in  terms  of  the  standard 
normal  distribution: 

7  =  P(fi  -  ka  <  fi—  Kpo) 


where  as  above 

6  =  K0 

The  problem  is  to  determine  a  function  MSf,S2)  so  that  (  33)  is  ap¬ 
proximately  satisfied  for  all  a f  and  <r\.  If  tolerance  limits  on  the  median  are 
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desired,  then  8  =  0  and  the  results  of  Welch  (1947)  and  Aspin  (1948)  may 
be  used  directly.  If  8  is  not  zero,  the  idea  behind  the  Welch- Aspin  derivation 
may  still  be  applied,  although  the  algebra  is  considerably  messier. 

The  Welch-Aspin  approach  makes  use  of  certain  differential  operators 
in  developing  an  ‘asymptotic’  series  for  h.  The  same  approach  may  be 
used  here,  but  for  first  order  calculations  the  additional  formalism  is  not 
justified  in  terms  of  algebraic  simplifications.  For  this  reason,  the  discussion 
below  consists  of  a  straightforward  Taylor  series  derivation.  Of  course,  both 
methods  must  give  the  same  answer,  and  this  has  been  used  to  provide  a 
check  on  the  calculations. 

Begin  by  rewriting  (  33)  as 


£(•!>(/,-, +  (7)1  =  7, 

(35) 

w'here 

=  ms;,s?)_6 

Oxisfn 

(36) 

Expand  h 

in  a  series  of  inverse  powers  of  n,, 

h  =  h.Q  -f  hi  +  0(l/m2), 

(37) 

where 


m  =  min(ni,  U2). 


Up  to  terms  of  second  order  in  h  we  have  that 


[  ho(SlS22)  K go 

Oil^/n  o\f\/n  cr\/y/n 


(38) 


where  we  have  substituted  (  34)  for  6. 

For  the  zeroth  order  approximation  we  approximate  h(S2,Sl)  by 


h0(Sl,Sj)  «  h0{a],a\) 


and  we  have,  U  =  0  and 


or 


Mfri'gf)  -  KqQ 

0\jy/n 


(39) 


ho(SlSl) 


A%5, 

v/H 


+  KgS. 


(40) 


9 


For  the  first  order  expression,  we  approximate  h  by 


ho(SlSi)  +  «  ho(SlSl)  +  filial)  = 


M1+(*  *)]  +  Js  [‘  +  (»,  "  ‘). 

j+At^i.^a)  (41) 

then 

U  =  K-,Ui  +  k0u2  + 

o  l/v/n 

(42) 

where 

°\ 

(43) 

and 

U*-<ril>/z(v  O’ 

(44) 

Let  X/  denote  a  random  variable  with  /  degrees  of  freedom  and  define 


(45) 


for  i  =  1,2.  The  [/;  can  be  expressed  in  terms  of  the  Vn,  as  follows: 


Ui  =  (1  + Vni)1/2-  1, 


u2  = 


<?1  IVn 


1/2 


(46) 

.(47) 


After  expanding  the  square  roots  in  (  46)  and  (  47)  in  power  series, 
one  can  readily  obtain  approximations  to  the  first  two  moments  of  the  U, 
suitabie  for  first  order  calculations: 


E(U  \) 

ss 

1 

4ni  ’ 

(48) 

E(U?) 

w 

1 

2ni  ’ 

(49) 

E(U2) 

« 

_  l  [fi  +  2l 
4  ("1  ri2. 

(50) 

e(u2) 

ss 

a  ffli  |  a2' 

2a\/y/n  [nj  TI2J  ’ 

(51) 
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and 


E -  /,/VC*!,- 

(52) 

where 

ofn1/2 
fll  ^  (T3  J2 

(53) 

and 

«,  =  "yi  /-’(v*. 

o  i  a 6 

(54) 

The  next  step  is  to  expand  the  normal  cdf  about  A%,  so 
be  replaced  by  the  following  approximation: 

that  (35)  may 

£[$(ff,  +  tf)]  =  7«*(Jir,) 

+<HKJE(U)  -  K^(KJE(U2)/ 2, 

(55) 

where  <£(■)  denotes  the  standard  normal  density.  The  expectation  of  U  can 
be  determined  immediately  from  (  42),  (  48)  amd  (  50).  Since 

E(l°)  =  K*E{Ul)  +  KlE{Ul)  +  2K0K\E(UlU2 )  +  0(l/m2),  (56) 


we  need  only  substitute  (  49),  (51)  and  (  52)  into  (  56)  in  order  to  complete 
the  evaluation  of  (  55). 

To  complete  these  calculations,  solve  (  55)  for  hi{o\,cr\)  (note  that  ht 
appears  through  E(U)),  replace  each  occurrence  of  <rf  or  o2  with  52  or 
S 2  respectively  (i  =  1,2)  and  divide  hi(S2,S2)  by  5  to  finally  obtain  the 
tolerance  limit  factor  k.  The  terms  of  k  may  then  be  rearranged  to  reveal 
their  structure.  The  following  expression  for  k  is  one  possibility: 


k 


r,  ,  K-yW  ,  W  \K^K*  +  1) 

K,  +  ~7r  +  47i[  ^ 

2 KfiKlVIW  ,  KpK-flW2 

+ - + - 

"i 

KpVlW3  KjK^IjJ  -  1  )2W2 

+  Til  +  «2  Q2 

,  K0{J  -  1  )2y/I\V3' 

+  n2Q2 


(57) 
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where 


«,s(i  +  (M)/g) 


(58) 


-1/2 


and 

g  s  S.  (59) 

The  coverage  probability  for  the  above  approximation  as  a  function  of 
the  population  variance  ratio  is  plotted  in  Figure  2  for  a  (.90,  .95)  tolerance 
limit  and  J  —  5.  Note  that  for  many  batches  the  series  solution  performs 
well,  though  for  few  batches  it  is  anticonservative. 


6  An  Alternative  Solution  for  Unknown  R 


For  small  samples,  the  first  order  approximation  developed  above  may  not 
be  adequate,  and  higher  order  calculations  are  clearly  prohibitive.  An  al¬ 
ternative  approach  is  to  view  the  problem  as  an  integral  equation,  following 
Trickett  and  Welch  (1954). 

If  one  defines 

r  =  JR  +  1  (60) 

then  (  24)  may  be  written  as 


Tni+n7  ^(r)(n!  +  +  =  7, 


(61) 


where 


and  B  is  defined  in  (  14).  The  parameter  r  may  be  estimated  by  the  mean 
square  ratio  (  59): 

Q  =  rFB1,nj,  (63) 


where  Fni  tnj  denotes  a  random  variable  having  an  F  distribution  with  nj  and 
«2  degrees  of  freedom.  The  tolerance  limit  problem  reduces  to  determining 
a  function  k(Q)  such  that 


12 


for  all  values  of  r  >  1,  where  Z,  Vj,  and  Y2  are  as  in  Section  3.  This  is 
equivalent  to  the  integral  equation 


VT(k)  =  Eu 


[m+nj  ^(g)(«i  +  n2)1/2y/^ 


1  _  y  i 

+  LJL,6(t) 


where 


Q  = 


rn2X 

ni(l-X) 


=  7> 

(65) 

(66) 


and  the  expectation  is  with  respect  to  a  beta  density  with  parameter  v  = 

(ni/2,n2/2). 

In  Section  5,  we  derived  an  approximation  to  k(Q)  which  we  label 
here  ko(Q).  We  intend  to  improve  this  approximation  by  replacing  it  by 
ki(Q)  =  k0(Q)  +  tp(Q)  and  we  will  employ 


k(c,Q)  =  h(Q)  +  £i/>(Q), 


(67) 


noting  that  k^(Q)  =  £(1,Q). 

Expanding  VT  [&(e,Q)]  in  a  Taylor  series  gives  the  first  order  approxima 
tion  dV 

7(0  =  Vr  [k0(Q)  +  erP(Q)]  *  Vr(k0(Q))  +  <  r 


de 


t=0 


(68) 


Employing  this  result,  the  equation  leading  to  our  next  approximation  may 
be  written  as 


7  =  7(1)  = 
VT(ko)  +  E„ 


'^ni+nj 


'P(Q)(ni  +  n2)1/2^7~i  +  ^"7^ 
MO)(ni  +  «2)1/2/ttt  + 


(69) 


where  denotes  the  noncentral  t  density.  The  noncentral  t  density 

with  /  degrees  of  freedom  and  noncentrality  parameter  6  may  be  calculated 
by  means  of  the  following  formula  (Odeh  and  Owen,  1980,  p.  272): 


T/(z,S) 


(70) 
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Since  there  are  computer  subroutines  available  for  determining  the  non- 
central  t  cdf  (see,  e.g.,  Griffiths  and  Hill,  1985),  (  70)  is  very  useful  for 
computation. 

The  first  term  on  the  right  hand  side  of  (  69),  VT(fco)>  may  be  evaluated 
numerically  for  given  r  since  ko(Q)  is  a  known  function. 

The  second  term  can  not  be  evaluated  without  knowing  ip.  To  simplify 
matters  we  shall  pretend  that  ip{Q)  assumes  a  constant  value  and  can  be 
factored  out  of  the  expectation.  The  Trickett- Welch  approach  consists  of 
replacing  ip(Q)  by  rp(qo)  where  qo  is  that  value  of  Q  corresponding  to  the 
mean  xo  of  the  beta  random  variable  X,  i.e. 


rn2x  o 

q°  ~  ni(l  -  x0) 


rn2ni/(ni  4-  rt2) 
ni  n2/(nj  +  n2) 


Thus  we  have 

7  *  K(k0)  +  V’(l")l/lr(^o) 

where 


Vir(feo)  = 

Ev 


(m  +  n2),/a^  ~  + 

•^ni+n2  ^o(<?)(ni  + 


(71) 

(72) 

(73) 


and 


xp(r) 


7  ~  Vrih) 

VlT(ko)  ’ 


ki  (r)  =  k0{r )  +  ip(r ). 


(74) 

(75) 


Our  numerical  integration  depends  on  the  use  of  fco  for  a  mesh  of  Q  or  r 
values  from  0  to  oo.  We  can  compute  fci  for  the  same  mesh.  This  fci(£?)  can 
be  used  as  the  fco(Q)  for  the  next  iteration. 

The  approximation  which  allows  t p(Q)  to  be  removed  from  the  integrand 
is  crude.  It  is  certainly  not  obvious  that  this  procedure  will  provide  any  im¬ 
provement  on  the  first  approximation.  In  fact,  if  one  were  to  implement  the 
algorithm  presented  in  this  section  on  a  computer,  one  would  see  that  the 
coverage  probability  improves  only  slightly  before  the  successive  approxi¬ 
mations  begin  to  diverge.  By  employing  this  method  as  improved  in  the 
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following  section,  we  can  do  much  better.  The  improvements  of  Section 
7  result  in  an  algorithm  which  provides  a  solution  to  the  tolerance  limit 
problem  of  unknown  R  that  is  (for  practical  purposes)  exact.  All  mention 
of  results  of  applying  the  Trickett- Welch  method  in  this  paper  refer  to  the 
modification  to  be  discussed  in  the  next  section. 


7  A  Modification  of  the  Trickett- Welch 
Approach 


In  order  to  obtain  useful  results  from  the  integral  equation  approach  of 
Section  6,  it  is  necessary  to  improve  the  rough  approximation  by  which 
the  unknown  function  0  is  removed  from  the  the  expectation  in  (  69).  The 
technique  by  which  this  approximation  is  improved  is  based  on  evaluating 
ip{Q )  for  a  value  q\  of  Q  corresponding  to  xx  of  X  where  the  integrand  of  V\r 
in  (  73)  achieves  its  maximum,  instead  of  xo,  the  mean  of  the  beta  density, 
which  may  be  close  to  where  this  density  achieves  its  maximum. 

For  any  value  of  r,  we  can  determine  the  desired  desired  peak  xx(r) 
numerically,  and  define 

_  rn2xi(r) 

91  “  nx(l  -  xi(t))‘ 


It  is  fortunate  that  in  our  tolerance  limit  problem,  the  value  Xi(r)  is  nearly 
independent  of  r.  Thus,  ip{Q)  can  be  evaluated  at  or  very  nearly  at  a 
specified  grid  of  qi  values  by  adjusting  r  after  the  nearly  constant  value  xj 
of  xj(r)  is  approximated  for  a  typical  r  value. 

One  difficulty  with  the  above  proposal  arises  from  the  fact  that,  strictly 
speaking,  r  should  only  be  taken  to  be  greater  than  one,  in  which  case  the 
range  of  q\  values  is  from  n2Xi/[nx(l  -  xj)]  to  oo  instead  of  from  0  to  oo  as 
is  required  for  the  numerical  integration.  Since  n2xj/[nx(l  -  xj)]  turns  out 
to  be  relatively  small,  we  translate  the  value  of  gx  by  this  amount,  so  that 
the  range  of  q  values  will  be  0  to  oo .  In  other  words,  we  replace  ko(q\ )  in 
the  approximation 


7  =  7(1)  =  Vr(^o(9i))+  0(?t)^r(fco(?i))  (77) 


by  k0{q\  -  n2ii/[ni(l  -  Xj)]}.  After  this  approximation  is  carried  out,  the 
method  can  be  iterated  using 


nx(l  -  xj) 


n2xj 

ni(l  _  xi) 


(78) 
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to  replace  ko. 

With  each  iteration  the  value  of  the  constant  ij  is  likely  to  change  and 
should  be  recalculated. 

The  above  simple  improvement  of  the  integral  mean  value  theorem  ap¬ 
proximation  underlying  the  Trickett-Welch  approach  enables  one  to  calcu¬ 
late  tolerance  limit  factors  which  provide  very  nearly  the  nominal  coverage 
probability  even  for  few  batches  and  small  batch  size. 

The  Trickett-Welch  approach,  possibly  with  modifications  similar  to 
those  discussed  in  this  section,  promises  to  be  applicable  to  other  problems, 
two  of  which  are  considered  below. 


8  Other  Applications  of  the  Trickett-Welch 
Approach 

The  Trickett-Welch  approach  can  be  applied  to  a  wide  range  of  problems  of 
inference  in  the  presence  of  a  nuisance  parameter.  Two  examples  of  such 
problems  are  outlined  in  this  section. 


8.1  Confidence  Intervals  for  the  Population  Mean 


For  the  random  effects  model  of  Section  1  a  two  sided  confidence  interval  for 
the  population  mean  is  desired  which  attains  nearly  the  nominal  confidence 
level  whatever  the  population  gintraclass  correlation 


P  = 


'k 


R 


+  <jg  n  +  r 


(79) 


Let  Z?i(  )  be  an  unspecified  function  of  the  mean  square  ratio  (  59).  With 
notation  as  in  Sections  5  and  6  we  have 


7  =  P({j.  -  D\S  <  m  <  A  +  D\ S).  (80) 

This  is  easily  shown  to  be  equivalent  to  the  integral  equation 

1  +  7  -  F 
—  ~E 

where  the  expectation  is  with  respect  to  the  distributions  of  the  mean 
squares.  The  methodology  presented  in  this  paper  can  be  applied  directly. 
An  important  feature  of  this  example  is  that  it  provides  a  method  which, 
for  a  particular  simple  situation,  avoids  the  problematic  issue  of  ‘when  to 
pool’. 


f  (Di(S?/Sl)S\ 


(81) 
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8.2  Testing  the  Equality  of  Two  Normal  Percentiles 

Another  thinly  disguised  version  of  the  Behrens-Fisher  problem  is  the  prob¬ 
lem  of  testing  the  equality  of  two  normal  percentiles  where  population  means 
and  variances  are  unknown.  Two  statistics  for  performing  such  a  test  are 
proposed  by  Cox  and  Jaber  (1985).  These  tests  require  simulation  in  order 
to  obtain  approximate  critical  values  for  the  test  statistics.  The  method 
outlined  below,  though  its  properties  have  yet  to  be  examined,  requires  no 
Monte-Carlo  tables. 

We  wish  to  test  equality  of  the  100/3t/i  percentiles  of  two  normal  popu¬ 
lations  on  the  basis  of  simple  random  samples  from  each  population.  That 
is,  we  are  interested  in  testing  the  null  hypothesis 

H0  :  +  Kpffi  =  H2  +  K 00 2  (82) 

against  the  alternative 

H\  :  /xj  -f  KpiTi  ^  H2  +  K 002 ,  (83) 

where  w  and  a  are  the  population  mean  and  standard  deviantions  for  i  = 

1.2  and  Ii'p  is  defined  in  (  7). 

Denote  the  sample  means  and  variances  by  X{  and  Sf  and  let  the  sample 
sizes  be  n,.  Define  the  statistic 

T  =  (Xi  +  I\0S\  )-(X2  +  KpS2)  (84) 


We  propose  to  reject  the  null  hypothesis  when  |T|  is  sufficiently  large.  A 
function  Z?2  which  provides  such  a  test  (of  size  1  -  7)  can  be  shown  to 
satisfy  the  following  integral  equation,  where  as  above  the  expectations  are 
with  respect  to  the  mean  squares: 


^  /  Dj  -  KpjSi  -  S2  -h  (?2  -  04 ) 

—  D2  —  Kp{S\  —  S2  +  <7  2  —  <?i)\ 
Ja\!m  +  a\ln2  ) 


which  is  in  a  form  to  which  the  techniques  of  this  paper  may  be  applied. 


9  The  Distributions  of  the  Tolerance  Limits 


Once  the  function  k  has  been  determined  it  is  straightforward  to  calculate 
the  cumulative  distribution  function  of  the  tolerance  limit.  It  is  obviously 
preferable  to  compare  distributions  of  confidence  bounds  rather  than  merely 
confidence  levels,  and  we  make  such  a  comparison  in  this  section. 

We  consider  a  (.90, .95)  lower  tolerance  limit  for  a  normal  population 
with  tenth  percentile  zero  and  variance  one.  In  Figure  3,  the  cumulative 
distributions  for  I  =  J  =  5  of  the  proposed  tolerance  limit  are  presented  for 
various  values  of  the  intraclass  correlation  p  =  R/(R  +  1). 

Note  that  all  of  the  curves  pass  very  nearly  through  (0,  .95),  indicating 
the  striking  success  that  we  have  had  at  removing  the  nusiance  parameter, 
even  for  as  few  as  five  batches.  As  the  intraclass  correlation  is  increased  the 
random  effects  sample  goes  from  behaving  essentially  like  a  single  sample  of 
size  a  =  IJ  when  p  =  0  to  being  equivalent  to  a  single  batch  of  cize  I  when 
P  =  1- 

In  Figure  4  three  cdfs  are  plotted,  corresponding  to  the  Mee-Owen 
method,  the  proposed  method  and  the  the  solution  for  known  R.  The 
intraclass  correlation  is  taken  to  equal  zero  and  the  sample  size  is  again 
7  =  7  =  5.  Note  that  the  result  based  on  the  Trickett- Welch  approach  is 
clearly  preferable  to  the  Mee-Owen  solution  and  doesn’t  fare  too  badly  when 
compared  to  the  known- R  solution. 

10  Discussion 

The  situation  of  primary  interest  to  the  aircraft  industry,  (.90,  .95)  lower  tol¬ 
erance  limits,  is  used  here  for  illustration.  Four  methods  have  been  presented 
in  this  paper:  the  Mee-Owen  method  (Section  2),  a  modified  Mee-Owen 
method  (Section  2),  a  method  based  on  the  Welch-Aspin  series  (Section  5), 
and  a  method  based  on  an  integral  equation  (Sections  6,  7).  The  coverage 
probability  functions  corresponding  to  these  methods  are  numbered  1-4  in 
Figure  5  for  five  batches  each  of  size  five. 

The  integral  equation  approach  virtually  removes  the  nuisance  parameter 
from  the  problem.  The  Mee-Owen  method  has  the  disadvantage  of  being 
substantially  conservative  when  the  variance  ratio  is  small. 

Only  a  slight  reduction  in  this  conservatism  has  resulted  from  the  mod¬ 
ification  of  the  confidence  level  of  the  variance  ratio  estimate  (Section  2, 
Table  2). 
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The  Welch- Aspin  series  solution  is  clearly  not  suitable  for  as  few  as  five 
batches,  as  discussed  in  Section  5,  However,  it  is  easy  to  compute  and 
provides  an  adequate  starting  function  for  the  iterative  integral  equation 
approximation  (  69). 

From  the  rescaled  plot  of  the  coverage  probability  function  for  the  inte¬ 
gral  equation  solution  (Figure  6)  it  can  be  seen  that  for  R  >  1  the  actual 
coverage  probability  differs  from  .95  by  no  more  than  ±.00005.  This  small 
difference  can  be  attributed  to  the  limited  accuracy  of  the  numerical  integra¬ 
tion.  For  R  <  1,  however,  the  difference  in  the  actual  and  nominal  coverage 
probability  increases  substantially,  but  never  does  it  reach  a  magnitude  that 
warrants  concern  for  applications. 

Figure  6  illustrates  the  convergence  of  the  Trickett-Welch  approach  for 
various  values  of  the  intraclass  correlation.  Note  that  for  practical  purposes 
ten  iterations  is  adequate,  although  some  slight  improvement  may  result 
from  considering  more  iterations. 

1 1  Example 

The  data  in  Table  3  are  a  pseudo-random  sample  of  25  from  a  normal  dis¬ 
tribution  with  mean  50  and  standard  deviation  10.  These  data  have  been 
arbitrarily  grouped  into  five  batches  of  five.  By  fitting  a  one-way  random 
effects  model  to  these  data  one  obtains  : 


MSb  =  89.88,  MSe  =  158.6, 

(86) 

fi  =  48.30,  b\  -  144.9. 

(87) 

A  (.90,. 95)  lower  tolerance  limit  is  of  the  form 

T  =  A  ~  K*x- 

(88) 

For  the  method  of  Mee  and  Owen  (1983)  K  =  1.90.  If  the  Mee-Owen 
method  is  modified  as  suggested  in  in  Section  2,  then  K  only  decreases 
to  1.89.  The  series  solution  of  Section  5  gives  K  -  1.78  and  the  integral 
equation  of  Section  6  results  in  K  =  1.83.  The  tolerance  limit  estimates  are, 
respectively,  25.42,  25.54,  26.82  and  26.29.  These  values  may  be  compared 
with  the  tolerance  limit  estimate  for  the  pooled  data,  which  is  26.00. 
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12  Conclusion 


One-sided  tolerance  limits  for  random  effects  models  is  a  topic  of  consid¬ 
erable  importance  in  engineering  statistics.  The  purpose  of  this  paper  has 
been  to  consider  this  tolerance  limit  problem  from  the  point  of  view  of  the 
Welch  interpretation  of  the  Behrens- Fisher  problem.  This  approach  leads  to 
a  method  which  provides  essentially  the  nominal  coverage  probability  what¬ 
ever  the  value  of  the  nuisance  parameter.  We  have  demonstrated  that  in 
addition  to  excellent  coverage  properties,  the  distribution  of  the  proposed 
tolerance  limit  compares  favorably  with  an  exi  s  t(n g,p roc ed ure  and  with  an 
exact  solution  for  known  nuisance  parameter.  /rT- 
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Coverage  probability  for  a  (.90,  .95)  lower  tolerance 
limit  using  the  pooled  data  method  as  a  function  of 
the  variance  ratio. 


VARIANCE  RATIO  (R) 


Coverage  probabilities  for  (.90,  .95)  tolerance 
llaits  using  the  first  order  Welch  aeries  method 
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Distributions  of  (.90,  .95)  lower  tolerance 
limits  obtained  by  the  Trickett-Welch  method 
for  various  values  of  the  intraclass  correlation. 
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Appendix 

FORTRAN  source  code  listings 

Listed  below  are  the  routines  used  to  perform  the  calculations 
in  this  paper.  All  of  the  required  software  is  listed  with 
the  following  exceptions: 

1)  Routines  in  the  IMSL  library 

2)  TEKTRONIX  PLOT-10  graphics  subroutines 

3)  The  noncentral-t  distribution  with  non-integer 

degrees  of  freedom  ('TNC',  Algorithm  AS  243, 

Applied  Statistics  (1989)  v.  38) 

4)  Routines  called  by  'TNC'  above,  all  of  which  are 

in  Griffiths  and  Hill  (1985) . 

The  routines  listed  in  this  appendix  are  available  in 
computer  readable  form  at  no  charge  from  the  author.  Send  a 
floppy  disk  (IBM-PC)  or  a  magnetic  tape  for  a  copy  of  the 
source  files. 

This  code  is  a  prototype  intended  as  a  research  tool. 

It  is  not  suitable  in  the  present  form  for  general  purpose 
use . 


Main  programs : 

TRICK  —  Program  to  calculate  tolerance  limit  factor 
by  the  modified  Trickett-Welch  approach  for 
a  balanced  nested  mixed  model,  a  simple 
generalization  of  the  model  considered  in 
this  paper. 

PLTCDF  —  Program  to  calculate  and  plot  distribution 

functions  for  tolerance  limits.  This  program 
uses  as  input  tolerance  limit  factor  files 
created  by  program  'TRICK' . 

COVRGE  —  Program  to  calculate  the  coverage  probability 
vs.  intraclass  correlation  functions  from 
tolerance  limit  factor  files  created  by 
program  'TRICK' 

Subroutines : 


EVCDF 

FNCK 

FNCN 

FNCR 

FNCS 

FNCY 

FNCZ 


Routine  to  evaluate  the  cdf  of  a  tolerance 
limit.  Called  by  ' TLMCDF' . 


Function 

called 

by 

root 

finder 

in 

' INVNCT' . 

Function 

called 

by 

root 

finder 

in 

' INVSPN' 

Function 

called 

by 

root 

finder 

in 

'  KR'  . 

Function  called  by  maximization  routine  in 
' FSUP' . 

Function  called  by  numerical  integration 
subroutine  in  'GENDS2' . 

Function  called  by  numerical  integration 
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subroutine  in  ' EVCDF' . 

GENDS2  —  Function  to  calculate  cdf  of  a  generalized 
noncentral-t  random  variable. 

INISPL  —  Subroutine  to  initialize  spline  interpolation 

of  tolerance  limit  factor. 

INIT  —  Initialization  routine  for  'TRICK'. 

INTEQN  —  Top-level  subroutine  for  solving  Trickett-welch 
integral  equation.  Called  by  'TRICK'. 

INVNCT  —  Subroutine  to  determine  noncentral-t  quantiles. 

INVSPN  —  Subroutine  to  perform  inverse  spline  interpolation. 

Called  by  'MESH' . 

KFACT  —  Subroutine  to  determine  tolerance  limit  factor  for 

a  simple  random  sample  from  a  normal  distribution. 

KMO  —  Subroutine  to  calculate  the  Mee-Owen  tolerance 

limit  factor.  Since  the  Satterthwaite  degrees 
of  freedom  need  not  be  an  integer,  it  is  because 
of  this  routine  that  the  Applied  Statistics 
subroutine  'TNC'  is  used.  For  integer  degrees  of 
freedom  the  IMSL  routine  is  adequate. 

KR  —  Routine  to  determine  the  tolerance  limit  factor 

for  known  variance  ratio. 

KSPLN  —  Spline  interpolation  for  tolerance  limit  factor. 

MESH  —  Subroutine  to  improve  the  mesh  of  nuisance 

parameter  values.  Initially,  the  Welch  series 
is  evaluated  for  equally  spaced  values  of  the 
ratio  of  mean  squares.  A  spline  is  fit  to  this 
function,  the  ordinate  is  divided  into  equal 
intervals,  and  the  spline  is  inverted  to  provided 
new  abscissa  values  which  will  be  closer  together 
where  the  function  has  a  larger  derivative.  This 
new  mesh  is  used  for  all  future  iterations. 

NCTD1N  —  Called  by  ' NCTDRV' .  Noncentral-t  density. 

NCTD2N  —  Called  by  'NCTDRV' . 

NCTD3N  —  Called  by  'NCTDRV' . 

NCTD4N  —  Called  by  'NCTDRV'. 

NCTD5N  —  Called  by  'NCTDRV'. 

NCTDRV  —  Subroutine  to  recursively  calculate  any  of  the 

first  five  derivatives  of  the  noncentral 
t  distrivution.  Derivatives  beyond  the  first 
are  not  used  in  the  paper. 

NEXTK  —  Subroutine  which  calculates  the  next  iteration 

of  the  modified  Trickett-Welch  procedure.  Called 
by  'INTEQN',  which  is  called  by  program  'TRICK'. 

SUP  —  Function  to  find  the  maximum  of  the  Trickett-Welch 
integrand  in  order  to  improve  integral  mean  value 
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approximation  as  discussed  in  Section  7.  Called  by 
'NEXTK' . 

TLMCDF  —  Subroutine  to  calculate  the  distribution  of  the 

tolerance  limit.  Called  by  'PLTCDF'  and  'COVRGE'. 

WELCH  —  Function  to  calculate  the  first  order  Welch  series 
approximation . 

XINT  —  Function  to  evaluate  the  two  integrands  for  the 
modified  Trickett-Welch  algorithm.  Called  by 
'NEXTK'  and  'FNCS' . 
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program  covrge 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 


c 

c 

c 

c 


10 

20 


30 

c 

c 

c 

c 
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Program  to  determine  points  on  coverage  probability  vs. 
intraclass  correlation  curves  for  tolerance  limit  factors. 
This  program  uses  as  input  files  created  by  'TRICK'. 

**  Note:  This  program  sends  extensive  output  to  the 

terminal  which  may  be  used  for  later  plotting. 

On  many  computers  this  output  will  have  to  be 
redirected  to  a  file. 

real  x(500),  xk(500),  xdx(100),  cov(100) 

integer  n(3) 

character*20  flenme,  file2 

character*3  citer 

character*l  ans 

common  /kw/  known,  rho,  xknown 

data  lfn  /10/ 

—  Filename  (specified  in  'TRICK') 
write  (*,*)  'filename  ?' 

read  (*,'(a20)')  flenme 

—  Number  of  points  at  which  coverage  probability  is  to 
be  determined 

write  (*,*)  'points  per  curve  ?' 
read  (*,*)  nrho 

—  Indices  for  tolerance  limit  factor  files.  These  files 
are  output  from  'TRICK'.  The  index  corresponds  to  the 
iteration,  it  follows  the  hyphen  in  the  file's  name, 
write  (*,*)  ' l=range  of  indices,  2=individual  indices  ?' 
read  (*,*)  iopt 
if  (iopt  .eq.  1)  then 

write  (*,*)  'min  index,  max  index  ?' 
read  (*,*)  imin  ,  imax 

ndp  =  imax  -imin  +1 
do  10  i=l,  ndp 

xdx(i)  =  imin  +(i-l) 
continue 
else 

ndp  *  0 
continue 

write  (*,*)  'index  (0  to  quit)  ?' 
read  (*,*)  idxl 
if  (idxl  .eq.  0)  go  to  30 
ndp  «■  ndp  +1 
xdx(ndp)  -  idxl 
go  to  20 
continue 
end  if 


—  Get  parameters  corresponding  to  tolerance  limit  factor 
from  files  refered  to  above. 


(nf ix=l 

for  cased  treated 

write  (*,*) 

'  nfix,  i,  j  ?' 

read  (*,*) 

n 

k  *  n (2) 

1  -  n  (3) 

write  (*,*) 

'  p  ?' 

read  (*,*) 

P 
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write  (*,*)  '  confidence  ?' 
read  (*,*)  g 
c 

c  —  Intraclass  correlation  is  not  regarded  as  known. 

known  =  .false, 
c 

c  —  Loop  over  tolerance  limit  facotr  input  files. 

do  40  i=l,  ndp 
c 

c  —  Build  the  file's  name  using  the  root  name  'flenme'  and 

c  the  index  number  corresponding  to  the  current  iteration, 

write  (unit=citer,  fmt=' (al, i2) ' )  '-',int  (xcU(i)) 
if  (citer  (2:2)  .eq.  '  ')  citer  (2:2)  =  'O' 
lstnbk  =20 
50  continue 

if  (flenme  (lstnbk : lstnbk)  .eq.  '  ')  then 
lstnbk  =  lstnbk  -1 
go  to  50 
end  if 

file2  =  flenme  (1: lstnbk)  //citer 
c 

c  —  Loop  over  points  on  each  curve 

dr  =  1/  float  (nrho  -1) 

c 

c  —  Header  for  curve  (you  may  want  to  direct  succeeding  output 
c  to  a  file  for  later  plotting) . 


write 

(*, 

*) 

write 

(*, 

*) 

r 

Number  of  fixed  factors 

'  ,n(l) 

write 

(*, 

*) 

f 

Number  of  random  batches 

' ,n (2) 

write 

(*, 

*) 

$ 

Batch  size 

' , n  (3) 

write 

<*, 

*) 

* 

Quantile 

'  .P 

write 

<*, 

*) 

r 

Confidence 

'  >g 

write 

<*, 

*) 

t 

Trickett-Welch  iteration 

' , int (xdx(i) ) 

write 

<*, 

*) 

do  60 

j-1 

r 

nrho-l 

rho 

= 

(j 

-1)  *dr 

r 

rho 

/ <1 .  -rho) 

c 

c  —  Varaince  components  and  mean  corresponding  to  ooleranoe 

c  limit  factor.  The  mean  is  taken  to  equal  the  percentile. 

s2w  =1/(1  +r) 
s2b  =  r  *s2w 

xmu  =  anorin  (p)  *sqrt  (s2w  +s2b) 
c 

c  —  Subroutine  to  evaluate  the  coverage  probability 

call  tlmcdf  (xmu,  s2b,  s2w,  n,  p,  0.,  cov(i),  file2) 
c 

c  —  Write  out  intrclass  correlation  and  coverage  probability. 

write  (*,*)  rho,  cov(i) 

60  continue 

40  continue 

o 

stop 

end 

program  pltcdf 

c  - 

parameter  (maxpts  =  500,  maxcrv  ■  50) 
c 

c  Mark  Vangel,  Oct.  1988 

c 

c  Program  to  calculate  and  plot  the  distribution  of  a  tolerance 

c  limit  for  a  random  effects  model. 
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c 


c 


c 

c 

c 


c 

c 


1 

c 

c 

c 


c 

c 


c 

c 


c 

c 


c 

c 


c 


logical  known 

character*20  flenme 
character*l  ans 

dimension  cdf  (maxpts) ,  quant  (maxpts),  ipoint  (maxcrv) , 
$  crossx(2),  crossy (2) ,  n(3) 


common 

/kw/ 

known,  rho,  xknown 

ipoint 

(1)  = 

0 

write 

(*,*) 

'  number  of  points  per  plot  ?' 

read 

(*,  *) 

nq 

nplot 

0 

—  Get  parameters  which  are  constant 

over  plot 

(' 

nfix' 

»  number  of  fixed  effects 

in  nested 

write 

(*,*) 

'  nfix,  i,  j  ?' 

read 

(*,*) 

n 

write 

<*,*) 

'  p  ?' 

read 

(*,*) 

P 

write 

<*,  *) 

'  confidence  ?' 

read 

(*,*) 

g 

nfix  = 

n  (1) 

k 

n  (2) 

1 

n  (3) 

—  Loop  over  all  k-factor  files. 

write 

(*,*) 

'  0=brief  output,  Incomplete  output 

read 

<*,*) 

ibrief 

continue 


—  Specify  a  filename  for  tolerance  limit  factor 
(output  from  program  TRICK) 

write  (*,*)  '  Filename  (''k"  if  rho  is  known)  ?' 

read  (*,'(a20)')  flenme 

if  (flenme  .eq.  '  ')  go  to  2 

write  (*,*)  '  Rho  ?' 

read  (*,*)  rho 

if  (rho  .eq.  1)  rho  **  rho  -l.e-6 


—  Intraclass  correlation  known 
if  (flenme  .eq.  'k')  then 
known  =  .true, 
call  init  (p,  g,  n) 
else 

known  =  .false, 
end  if 


—  Var.  ratio,  var.  components 
r  *  rho  / (1  -rho) 

s2w  -  1.  /(I  +r) 
s2b  «  r  *s2w 

xmu  *  anorin  (p)  *sqrt  (s2w  +s2b) 

—  Pointer  to  next  plot 
nplot  -  nplot  +1 

ipoint  (nplot+1)  “  ipoint  (nplot)  +nq 

—  Calculate  cdf  at  equally  spaced  points  within  specified  range 
write  (*,*)  '  Range  of  values  for  cdf  ?' 

read  (*,*)  qmin,  qmax 
dq  “  (qmax-qmin) / (nq-1 . ) 


—  Header  for  plot 


percentile 


write 

<*, 

*) 

/ 

Distribution  of  ',100*(l-p),  '  percentil 

write 

(*, 

*) 

r 

from  ' ,  qmin,  '  to  ' ,  qmax 

write 

<*, 

*) 

r 

Mean  —  ',xmu,  '  Var.  components  -  ' , s2b. 

write 

<*, 

*) 

i 

Tolerance  limit  factor  file  =  ', flenme 

write 

(*, 

*) 

t 

Number  of  groups,  group  size  -  ' ,  k,  1 

write 

(*, 

*) 

9 

Confidence  level  =  ',  g 

write 

(*, 

*) 

c  —  Calculate  points  on  cdf 

do  20  i=l,  nq 

idx  =  (nplot  -l)*nq  +i 

quant  (idx)  =  (i-1)  *dq  +qmi.n 

call  tlmcdf  (xmu, s2b, s2w, n, p, quant (idx) , cdf (idx) , flenme) 
if  (ibrief  .eq.  1)  then 

write  (*,*)  quant  (idx),  cdf  (idx) 
end  if 

20  continue 

go  to  1 
c 

c  —  Now  plot  the  results 

2  continue 

write  (*,*)  'Plots  ?' 
read  (*, '  (al) ' )  ans 
if  (ans  .eq.  'y')  then 

write  (*,*)  '  Min  and  max  for  abscissa  ?' 
read  (*,*)  qmin,  qmax 

write  (*,*)  '  Min  and  max  for  ordinate  [0,0  for  default]  ?' 
read  (*,*)  omin,  omax 

c 

c  —  Initialize  PLOT-10  graphics 

call  initt  (960) 
call  binitt 

c 

c  —  Set  coordinate  ranges 

call  comset  (ibasex(ll),  qmin) 
call  comset  (ibasex(12),  qmax) 
if  (omax  .ne.  0.)  then 

call  comset  (ibasey(ll),  omin) 
call  comset  (ibasey(12),  omax) 
end  if 
c 

c  —  Produce  the  plots 

call  npts  (nplot  *nq) 
call  check  (quant,  cdf) 
call  npts  (nq) 
call  dsplay  (quant,  cdf) 
do  30  i=l,  nplot-1 

call  cplot  (quant  (i*nq+l),  cdf  (i*nq+l) ) 

30  continue 

c 

c  —  Plot  crosshairs  at  quantile  and  nominal  coverage  probability 

crossx  (1)  =■  xmu  -anorin  (p)  *sqrt  (s2w  +s2b) 
crossx  (2)  *=  crossx  (1) 
crossy  (1)  **  0. 
crossy  (2)  =  1. 
call  npts  (2) 

call  cplot  (crossx,  crossy) 
crossx  (1)  -  qmin 
crossx  (2)  “  qmax 
crossy  (1)  =  g 
crossy  (2)  =  g 
call  cplot  (crossx,  crossy) 
c 

c  --  Hardcopy  option 
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c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 


call  scursr  (ans,  il,  i2) 
if  (ans  .eq.  'p')  then 
call  hdcopy 
end  if 

call  finitt  (-1) 
go  to  2 
end  if 

stop 

end 

program  trick 
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Program  to  calculate  one  sided  tolerance  limits  by  the 
modified  Trickett-Welch  method. 

logical  restrt 
dimension  n(3) 
character*20  flenme,  rsfile 
common  /crs/  restrt,  rsfile 

--  Restart  capability  allows  restarting  from  a  previously 
computed  k  function, 
write  (*,*)  'Restart  (1  or  0)  ?' 
read  (*,*)  ires 
if  (ires  .eq.  1)  then 
restrt  =  .true, 
write  (*,*)  'Restart  file  ?' 
read  <*,'(a20)')  rsfile 
else 

restrt  =  .false, 
end  if 


—  Problem  parameters 

write  (*,*)  'Number  of  steps  for  ratio  of  mean  squares  ?' 
read  <*, *)  nstp 

write  (*,*)  'Filename  for  k-function  files  ?' 
read  (*,'(a20)')  flenme 
write  (*,*)  'Number  of  iterations  ?' 
read  (*, *)  niter 

write  (*,*)  '(1-quantile)  for  tolerance  limit  ?' 
read  (*,*)  p 

write  (*,*)  'Confidence  coefficient  for  lower  limit  ?' 
read  (*, *)  g 

write  (*,*)  'Number  of  fixed  effects  (for  nested  model)  ?' 
read  (*,*)  n  (1) 

write  (*,*)  'Number  of  random  batches  ?' 

read  (*,*)  n  (2) 

write  (*,*)  'Batch  size  ?' 

read  (*,*)  n  (3) 

call  inteqn  (p,  g,  n,  nstp,  flenme,  niter,  1.) 

stop 

end 
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c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 


c 

c 

c 


c 


c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 


subroutine  evcdf  (cum,  idfla,  idf2a,  cla,  c2a,  xncpa,  etaa,  serr) 
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Routine  called  by  ' TLMCDF' . 

double  precision  aerr,  error,  xl,  xh,  result,  fncz 
external  fncz 

—  Parameters  for  FNCZ 

common  /bl/  idfl,  idf2,  cl,  c2,  xncp,  con,  eta 
data  hf  /  .  5/ 

—  Double  precision  error 
aerr  -  serr 

—  Integration  rule 
irule  =  2 

—  Put  stuff  in  common 
idfl  =  idfla 

idf2  =  idf2a 
eta  =  etaa 
cl  =  cla 
c2  =  c2a 
xncp  =  xncpa 

con  =  alngam(hf * (idf l+idf2) )  -aingam(hf *idf 1)  -alngam(hf *idf2) 

—  Limits  of  integration.  Avoid  zero  and  one. 

xl  =  l.d-10 

xh  =  l.dO  -  l.d-10 

Do  the  integration. 

call  dqdag  (fncz,  xl,  xh,  aerr,  O.dO,  irule,  result,  error) 
cum  =  result 

return 

end 

real  function  fnck  (x) 


real*8  tnc 

common  /kcom/  xncp,  g,  idf,  df 

Called  by  root  finder  in  'INVNCT'. 

—  Noncentral  t  with  non-integer  degrees  of  freedom 
(Satterthwaite  d.f.  need  not  be  an  integer) 
fnck  «  tnc  (dble(x),  dble(df),  dble (xncp) ,  ifault)  -g 
return 
end 

real  function  fncn  (x) 


Called  by  root  finder  in  ' INVSPN' . 

common  /sp2/  y 
call  kspln  (x,  yl) 
fncn  -  y  -yl 
return 
end 

real  function  fncr  (x) 
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Mark  Vangel,  June  1986 


c 
c 

c  Routine  called  by  root  finder  in  subroutine  '  KR' . 

c 

common  /krl/  cl,  c2,  idfl,  idf2,  xkp,  xkg,  p,  g,  xncp 
fncr  -  g  -gends2  (x,  idfl,  idf2,  cl,  c2,  xncp,  l.e-7) 
return 
end 

real  function  fncs  (x! 

c  - 

c 

c  Called  by  maximization  routine  in  'FSUP' . 
c 

double  precision  xint 

common  /ca/  nfix,  k,  1,  eta,  xkp,  xkg,  idrv,  con 
c 

idrv  =  1 

fncs  =  -xint  (dble  (x) ) 

return 

end 

double  precision  function  fncy  (f) 

c  - 

c 

c  Mark  Vangel,  June  1986 

c 

c  Called  by  numerical  integration  subroutine  in  'GENDS2' . 

c 

implicit  real  (a-h,  o-z) 
double  precision  f 

common  /bl/  idfl,  idf2,  tval,  cl,  c2,  xncp,  con 
data  hf  /0.5/ 

fncy  -  dble  < (hf*idf2-l)  *log  (f)  +(hf*idfl-l)  *log  (1  -f ) ) 
arg  -  tval  *sqrt  (cl*(l.  -f)  +c2*f) 
if  (arg  .gt.  I.elO)  then 
tprob  -  1 . 

else  if  (arg  .It.  -I.elO)  then 
tprob  =  0 . 

else 

tprob  =  tndf  (arg,  idfl+idf2,  xncp) 
end  if 

fncy  =  dble (exp (fncy  +con)  *tprob) 

return 

end 

double  precision  function  fncz  (f) 

c  - 

c 

c  Called  by  ' DQDAG'  in  ' EVCDF' . 
c 

implicit  real  (a-h,  o-z) 
double  precision  f 

common  /bl/  idfl,  idf2,  cl,  c2,  xncp,  con,  eta 
data  hf,  one,  zero  /.5,  1.,  0./ 
c 

fncz  =  (hf *idf 2-one)  *dlog  (one  -f)  + (hf *idf 1-one)  *dlog  (f) 
x  **  eta  *idf2  *f  /((idfl  *(one  -f))) 
c 

c  —  Spline  interpolation  of  tolerance  limit  factor 

call  kspln  (x,  xk) 

arg  -  xk  *sqrt  (cl  *f  +c2  *(one  -f)) 
if  (arg  .gt.  I.elO)  then 
tprob  «  one 

else  if  (arg  .It.  -I.elO)  then 
tprob  -  zero 

else 
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o  o 


tprob  -  tndf  (arg,  idfl+idf2,  xncp) 
end  if 

fncz  -  dble  (exp  (fncz  +con)  *tprob) 
c 

return 
end 

real  function  gends2  (tvala,  idfla,  idf2a,  cla, 

$  c2a,  xncpa,  aerr) 

c  - 

c 

c  Mark  Vangel,  June,  1986 
c 

c  Evaluate  generalized  non-central  t  using  integral 

c  representation, 

c 

c  tvala 

c  idfla,  idf2a 

c  cla,  c2a 

c  xncpa 

c  serr 

c 

implicit  real  (a-h,  o-z) 

double  precision  aerr,  error,  xl,  xh,  result,  fncy,  rerr 
external  fncy 

—  Constants  for  fncy 

common  /bl/  idfl,  idf2,  tval,  cl,  c2,  xncp,  con 
c 

c  —  Constants  for  common  block 

idfl  =  idfla 
idf2  =  idf2a 
tval  =  tvala 
Cl  =  cla 
c2  -  c2a 
xncp  =  xncpa 
c 

c  —  Constants  for  numerical  integration 

aerr  =  serr 
rerr  =  O.dO 
hf  =0.5 
xl  =  l.d-10 
xh  =  l.dO  -  xl 

con  =  alngam(hf * (idf l+idf2) )  -alngam(hf*idf 1)  -alngam(hf *idf2) 

call  dqdags  (fncy,  xl,  xh,  aerr,  rerr,  result,  error) 

gends2  =  sngl  (result) 

return 

end 

subroutine  inispl  (lfn) 

c  - 

c 

c  Mark  Vangel,  Oct.  1988 

c 

c  Routine  to  initialize  a  spline  fit  to  data  read  from 

c  a  file.  The  unit  number  of  the  file  is  'lfn'.  The  tolerance 

c  limit  files  are  output  from  program  'TRICK', 

c 

common  /cb/  xs  (500),  xks  (500),  break  (500),  c  (4,  500),  nx 
c 

read  (lfn,  *)  nx 
do  10  i=l,  nx 

read  (lfn,  *)  xs  (i),  xks  (i) 

10  continue 

c 

call  csint  (nx,  xs,  xks,  break,  c) 


—  Argument  of  gen  net 

—  Degrees  of  freedom  for  chisquares 

—  Corresponding  coefficients 

—  Noncentrality  parameter 

—  Absolute  error  for  num.  integration 


All 


c 


c 

c 

c 

c 

c 

c 


c 


c 

c 


c 

c 


c 

c 


c 

c 

c 


c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


return 

end 

subroutine  init  (p,  g,  n) 


Mark  Vangel,  Oct.  1988 

Initialize  a  few  constants  that  will  be  needed  later. 

logical  known,  meowen 
dimension  n  (3) 

common  /ca/  nfix,  k,  1,  eta,  xkp,  xkg,  idrv,  con 
common  /kw/  known,  rho,  xknown 
common  /cc/  xkzero,  xkinf 
common  /cd/  diag 

data  one  /l./,  half  /0.5/,  rhoeps  /l.e-4/ 

nfix  *>  n  (1) 
k  =  n  (2) 

1  =  n  (3) 

—  Normal  quantiles, 
xkp  -  anorin  (p) 
xkg  =  anorin  (g) 

—  Log  beta  function, 

df 1  =  nfix  * (k  -one) 

df2  =  nfix  *(k  *(1  -one)) 
con  =  alngam  (half  * (df 1  +df2) )  - 
$  alngam  (half  *dfl)  -alngam  (half  *df2) 

—  k  values  for  zero  and  infinity 
call  kr  (n,  p,  g,  1.,  xkzero) 
call  kfact  (p,  g,  k-1,  xkinf) 

if  (diag  .eq.  one)  write  (*,*)  'init  :  xkzero,  xkinf  ', 

$  xkzero,  xkinf 

—  If  rho  is  known,  calculate  constant  k  value;  using  limit  for 
rho  =  one 

if  (known  .and.  abs  (rho  -one)  .gt.  rhoeps)  then 
t  =1  *rho  / (one  -rho)  +one 
call  kr  (n,  p,  g,  t,  xknown) 
else  if  (known)  then 
xknown  «  xkinf 
end  if 

return 

end 

subroutine  inteqn  (pa,  ga,  n,  nstpa,  flenme,  niter,  diaga) 


Mark  Vangel,  Oct.  1988 

Top-level  subroutine  to  compute  tolerance  limits  by  a 
modified  Trickett-Welch  procedure . 


P 

g 

n  (1) 
n  (2) 
n  (3) 
nstp 
flenme 


Probability  associated  with  quantile 
Confidence  associated  with  tolerance  limit 
Number  of  fixed  factor  levels  (nfix) 

Number  of  batches  (k) 

Batch  size  (1) 

Number  of  steps  for  nuisance  parameter 
Filename  for  output  of  tolerance  limit  factor 
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c 

c 

c 

c 


c 


c 


c 

c 

c 


c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 


c 

c 


20 

c 

c 


estimates 

niter  —  Number  of  iterations 
diag  —  =1  for  debug  lines  to  print 

parameter  (maxpts  =  500) 
logical  known,  restrt 

character*20  flenme,  rsfile 
character*30  file2 
character*3  citer 

dimension  xk  (maxpts),  xkl  (maxpts),  cvrate  (maxpts), 

$  xkstep  (maxpts) , 

$  x  (maxpts)  ,  n (3) 

common  /crs/  restrt,  rsfile 
common  /kw/  known,  rho,  xknown 

common  /ca/  nfix,  k,  1,  eta,  xkp,  xkg,  idrv,  con 
common  /cc/  xkzero,  xkinf 
common  /cd/  diag 

data  zero  /0./,  one  /l./,  two  /2./  ,  lfn  /10/,  rh  /20./ 

—  Initialization 
known  =  .false, 
nfix  =  n  (1) 

k  =  n  (2) 

1  =  n  (3) 

P  =  pa 
g  =  ga 
nstp  *  nstpa 
diag  =  diaga 

—  Degrees  of  freedom 
dfl  =  nfix  * (k  -one) 

df2  =  nfix  * (k  *(1  -one)) 

—  Constants  for  repeated  use 
call  init  (p,  g,  n) 

—  Initial  step  size  for  x 

dx  -  (l*rh  tone)  /(nstp  -one) 

iter  =  0 

—  First  guess  at  value  at  zero.  Note  that  Welch  result 
blows  up  at  zero,  hence  it  can't  be  used  here. 

x  (1)  =  zero 

xk  (1)  =  xkzero 

cvrate  (1)  =  g 

—  Values  at  infinity  :  exact 

x  (n3tp  +1)  -  one 

xk  (nstp  +1)  -  xkinf 

cvrate  (nstp  +1)  *=  g 

—  Option  to  continue  a  previous  calculation 
if  (restrt)  then 

open  (unit-10,  file-rsfile) 
do  20  i-1,  nstp  +1 

read  (10,  *)  x  (i),  xk  (i) 
continue 
else 

—  First  pass  at  Welch  series  :  equally  spaced  abcissas 
do  21  i-2,  nstp 
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x  (i)  «=  (i  -one)  *dx 

xk  (i)  -  welch  (x  (i) ,  xkg,  xkp,  n) 

21  continue 
c 

c  —  Write  out  first  pass  to  a  scratch  file  for  'mesh' 

open  (unit=lfn,  status**' scratch'  ) 
write  <lfn,  *)  nstp 
do  22  i=l,  nstp  +1 

write  (lfn,  *)  x  (i)  ,  xk  (i) 

22  continue 
c 

c  —  Find  mesh  which  gives  equal  spacing  in  y 

rewind  (lfn) 

call  inispl  (lfn) 
call  mesh  (x) 
close  (lfn) 

c 

c  —  Second  pass  at  Welch  series  :  equally  spaced  ordinates 

do  23  i=2,  nstp 

xk  (i)  *=  welch  (x  (i)  ,  xkg,  xkp,  n) 

23  continue 
c 

end  if 
c 

c  —  Loop  over  iterations 

do  30  i=l,  niter 
c 

c  —  Filename  for  output.  Iteration  number  appended  to  name, 

write  (unit=citer,  fmt=' (al, i2) ' )  ' , i 
if  (citer  (2:2)  .eq.  '  ')  citer  (2:2)  =  '0' 
lstnbk  =  20 
31  continue 

if  (flenme  (lstnbk : lstnbk)  .eq.  '  ')  then 
lstnbk  =  lstnbk  -1 
go  to  31 
end  if 

file2  *  flenme  (1: lstnbk)  //citer 
c 

c  —  Write  out  latest  results  to  file 

open  (unit=lfn,  file=*file2,  status®' new' ) 
write  (lfn,  *)  nstp 
do  40  j-1,  nstp  +1 

write  (lfn,  *)  x  (j),  xk  (j) 

40  continue 

c 

c  —  Initialize  the  spline  with  latest  results 

rewind  (lfn) 
call  inispl  (lfn) 
c 

c  —  Use  Trickett-Welch  to  get  improved  approximation  to  xk 

call  nextk  (n,  nstp,  x,  xk,  xkl,  xkstep,  cvrate,  p,  g) 
c 

c  —  Rewrite  the  current  approximation  with  coverage  rates 

rewind  (lfn) 
write  (lfn,  *)  nstp 
do  50  j“l,  nstp  +1 

write  (lfn,  *)  x  (j),  xk  (j),  cvrate  (j) 

50  continue 

c 

c  —  Update  current  approximation 

do  60  j-1,  nstp 

xk  (j)  -  xkl  (j) 

60  continue 

c 
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c  —  Now  do  it  all  over  again  . . . 

iter  =  iter  +1 
30  continue 
c 

return 

end 

subroutine  invnct  (ga,  dfa,  xncpa,  xl,  xh,  t) 

c  - 

c 

c  Mark  Vangel,  Dec.  1988 

c 

c  Subroutine  to  invert  the  noncentral  t  distribution,  the 

c  limits  'xl'  and  'xh'  contain  the  root  and  are  input  parameters, 

c 

common  /kcom/  xncp,  g,  idf,  df 
external  fnck 

data  aerr  /l.e-5/,  rerr  /l.e-5/ 
c 

g  =  ga 
xncp  =  xncpa 
df  =  dfa 
idf  =  df 
a  =  xl 
t  =  xh 

maxfn  =  250 

call  zbren  (fnck,  aerr,  rerr,  a,  t,  maxfn) 

return 

end 

subroutine  invspn  (xla,  xha,  ya,  x) 

c  - 

c 

c  Mark  Vangel,  Dec.  1988 

c 

c  Subroutine  'INVSPN'  performs  inverse  spline  interpolation, 

c  This  routine  is  called  by  'MESH" 

c 

common  /sp2/  y 
external  fncn 

data  zero  /0.0/,  eps  /l.e-5/ 
c 

xl  ~  xla 
xh  =  xha 
y  *  ya 
C 

errrel  =  eps 
errabs  =  zero 
maxfn  =  100 

call  zbren  (fncn,  errabs,  errrel,  xl,  xh,  maxfn) 

x  =  xh 

return 

end 

subroutine  kfact  (p,  g,  idf,  xk) 

c  - 

c 

c  Mark  Vangel,  June  1986 

c 

c  Subroutine  to  compute  tolerance  limit  factor  for  a 

c  simple  normal  sample, 

c 

data  one  /l./,  uplim  /25./ 
c 

xncp  -  anorin  (p)  *3qrt  (idf  +one) 

call  invnct  (g,  float (idf),  xncp,  one,  uplim,  t) 

xk  «  t  /sqrt  (idf  +one) 
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c 

return 

end 

subroutine  kmo  (i,  j,  p,  g,  xmsr,  fp,  xk) 

c  - 

c 

c  Mark  Vangel,  Dec.  1988 

c 

c  Calculate  the  Mee-Owen  tolerance  limit  factor, 

c 

data  zero  /0/,  half  / .  5 / ,  one  /1.0/,  xh  /25/ 
c 

c  —  upper  confidence  bound  on  variance  ratio 

dfl  =  i  -one 

df2  =  i  Mj-one) 

fconf  *  fin  (fp,  df2,  dfl) 

r  =  (xmsr*fconf  -one)  /real  (j) 

c 

c  —  noncentrality  parameter 

xb  =  (r  +one)  /  (j*r  tone) 

xncp  =  sqrt  (i*j  *xb)  *anorin(p) 
c 

c  —  Satterthwaite  degrees  of  freedom 

sdf  =  (r  +one) **2  / 

&  ( (r+one/ j) **2/ (i-1)  + (one-one/ j) /  (i*j) ) 

c 

c  —  noncentral  t  quantile 

call  invnct  (g,  sdf,  xncp,  one,  xh,  xk) 

c 

c  —  tolerance  limit  factor 
xk  =  xk  /sqrt  (i*j  *xb) 
c 

return 

end 

subroutine  kr  (n,  pa,  ga,  teta,  xk) 

c  - 

c 

c  Mark  Vangel,  June  1986 

c 

c  Routine  to  determine  tolerance  limit  factors  for  known 

c  variance  ratio  r. 
c 

c  n  —  number  fixed  effects,  batches,  batch  size 

c  pa  —  quantile 

c  ga  —  confidence  for  lower  tolerance  limit 

c  teta  —  eta=j*r+l  (known)  nuisance  parameter 

c  xk  —  returned  tolerance  limit  factor 

c 

dimension  n  (3) 
external  fncr 

common  /krl/  cl,  c2,  idfl,  idf2,  xkp,  xkg,  p,  g,  xncp 
c 

c  —  Constants  for  common  block 

eta  «  teta 
P  “  Pa 

g  -  ga 

xkg  -  anorin  (g) 

xkp  -  anorin  (p) 

if  (eta  .eq.  0.)  eta  -  .05 
c 

c  —  Degrees  of  freedom 

nfix  -  n  (1) 

i  -  n  (2) 

j  -  n  (3) 
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idfl  =  nfix  * ( i— 1 > 
idf2  -  nfix  * (i*  ( j-1) ) 
c 

c  —  Coefficients  for  linear  comb,  of  chi-squares 

cl  =  i*j  -1. 
c2  =  cl  /eta 
cl  =  cl  *i/(i-l.) 
c 

c  —  Noncentrality  parameter 

xncp  =  xkp  *sqrt  (i*(l.  + ( j-1 . ) /eta) ) 
c 

c  —  Root  finder 

aerr  =  l.e-5 
rerr  =  l.e-5 
a  =1. 

b  =10.  *xncp 

maxfn  =  100 

call  zbren  (fncr,  aerr,  rerr,  a,  b,  maxfn) 

xk  =  b 

return 

end 

subroutine  kspln  (eta,  xk) 

c  - 

c 

c  Mark  Vangel,  Oct.  1988 

c 

c  Spline  interpolation  of  tolerance  limit  factor, 

c 

logical  known 
dimension  m(3) 

common  /ca/  nfix,  k,  1,  deta,  xkp,  xkg,  idrv,  con 
common  /cb/  xs  (500),  xks  (500),  break  (500),  c  (4,  500),  nx 
common  /kw/  known,  rho,  xknown 
data  iord  /10/,  one  /1.0/ 
c 

c  —  Use  constant  value  if  rho  is  known 

if  (known  .and.  rho  .ge.  0)  then 
xk  =  xknown 
c 

c  —  Truncate  function  at  upper  limit  calculated 

else  if  (eta  .gt.  xs (nx) )  then 
xk  =  xks  (nx) 

else 

c 

c  —  Spline  interpolation 

xk  =  csval  (eta,  nx-1,  break,  c) 
end  if 
return 
end 

subroutine  me3h  (xnew) 

c  - 

common  /cb/  xs  (500),  xks  (500),  break  (500),  c  (4,  500),  nx 
c 

c  Mark  Vangel,  Dec.  1988 

c 

c  The  initial  mesh  of  abscissa  values  is  taken  to  be  equally 

c  spaced.  A  better  approximation  can  be  obtained  if  more  points 
c  are  taken  where  the  function  being  estimated  changes  most 
c  rapidly,  however.  Subroutine  'MESH'  takes  as  input  the 
c  initial  equally  spaced  mesh  and  the  Welch  approximation  at 
c  these  mesh  points.  The  ordinate  is  equally  divided  into 

c  intervals  and  the  Welch  approximation  provides  (via  inverse 

c  interpolation  in  ' INVSPN' )  the  corresponding  new  mesh  of 

c  abscissa  values.  This  new  mesh  is  used  for  all  succeeding 
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c 

c 


iteration . 


dimension  xnew  (1) 
data  one  /!./ 
c 

dlty  *  (xks (nx)  -xks(l))  / (nx-one) 
do  10  i-2,  nx-1 

y  *  ( i— 1 )  *dlty  +xks{l) 

call  invspn  (xs(l),  xs(nx),  y,  xnew(i)) 

10  continue 

xnew(l)  =  xs (1) 
xnew(nx)  -  xs (nx) 
return 
end 

subroutine  nctdln  (idf,  tval,  xncp,  densty) 

c  - 

c 

pi  -  tndf  (sqrt ( (idf +2 . ) /idf )  *tval,  idf +2,  xncp) 
p2  =  tndf  (tval,  idf,  xncp) 

c 

densty  =  (idf/tval)  * (pi  -p2) 

return 

end 

subroutine  nctd2n  (idf,  tval,  xncp,  drv) 

c  - 

c 

c  -  sqrt  ( (idf+2 . ) /idf ) 

call  nctdln  (idf+2,  c*tval,  xncp,  pi) 

call  nctdln  (idf  .  tval,  xncp,  p2) 

drv  =  (idf/tval)  * (pi  -p2) 

return 

end 

subroutine  nctd3n  (idf,  tval,  xncp,  drv) 

c  - - 

c 

c  -  sqrt  ( (idf+2 .) /idf ) 

call  nctd2n  (idf+2,  c*tval,  xncp,  pi) 

call  nctd2n  (idf,  tval,  xncp,  p2) 

drv  *=  (idf/tval)  *  (pi  -p2) 

return 

end 

subroutine  nctd4n  (idf,  tval,  xncp,  drv) 

c  - 

c 

c  =  sqrt  ( (idf+2 .) /idf ) 

call  nctd3n  (idf+2,  c*tval,  xncp,  pi) 

call  nctd3n  (idf,  tval,  xncp,  p2) 

drv  -  (idf/tval)  * (pi  -p2) 

return 

end 

subroutine  nctd5n  (idf,  tval,  xncp,  drv) 

c  - 

c 

c  -  sqrt  ( (idf+2 .) /idf ) 

call  nctd4n  (idf+2,  c*tval,  xncp,  pi) 

call  nctd4n  (idf,  tval,  xncp,  p2) 

drv  «■  (idf/tval)  *  (pi  -p2) 

return 

end 

subroutine  nctdrv  (k,  idf,  tval,  xncp,  drv) 

c  - 

c 

c  Mark  Vangel,  May  1986 

c 
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c  Evaluate  either  the  noncentral  t  cumulative  (k»0)  or 

c  the  kth  derivative  of  the  cumulative  with  respect  to  the 

c  argument  (k=l, 2, 3, 4, 5) .  derivatives  are  calculated  exactly 

c  in  terms  of  the  cumulative  by  means  of  a  recursion  formula, 

c 

if  (k  .eq.  0)  then 

drv  =  tndf  (tval,  idf,  xncp) 
else  if  (k  .eq.  1)  then 

call  nctdln  (idf,  tval,  xncp,  drv) 
else  if  (k  .eq.  2\  then 

call  nctd2n  (idf,  tval,  xncp,  drv) 
else  if  (k  .eq.  3)  then 

call  nctd3n  (idf,  tval,  xncp,  drv) 
else  if  (k  .eq.  4)  then 

call  nctd4n  (idf,  tval,  xncp,  drv) 
else  if  (k  .eq.  5)  then 

call  nctd5n  (idf,  tval,  xncp,  drv) 
end  if 
return 
end 

subroutine  nextk 

$  (n,  nstp,  x,  xkO,  xkl,  xkstep,  cvrate,  pa,  ga) 

c  - 

c 

c  Mark  Vangel,  Oct.  1988 

c 

c  Given  an  input  tolerance  limit  factor  xkO  and  the  parameters 

c  of  the  problem,  this  subroutine  calculates  the  next  approximation 

c  xkl  by  a  modified  Trickett-Welch  procedure, 

c 

c  n  (1) 

c  n  (2) 

c  n  (3) 

c  nstp 

c  x 

c  xkO 

c  xkl 

c  cvrate  — 

c  p  — 

c  g  — 

c 

dimension  x(l),  xkO  (1),  xkl  (1),  xkstep  (1),  cvrate  (1),  n(3) 
double  precision  tl,  th,  daerr,  drerr,  xiO,  xil,  errest 
common  /ca/  nfix,  k,  1,  eta,  xkp,  xkg,  idrv,  con 
common  /cd/  diag 
common  /sblk/  oldsup 
external  xint 

data  tl  /l.d-5/,  th  /l.dO/,  aerr/0.0/,  rerr/l.e-4/,  one  /l./, 

$  zero  /0./,  daerr  /l.d-5/,  drerr  /O.dO/ 
c 

c  —  Initialize  some  constants 
th  -  one  -tl 
nfix  «  n  (1) 

k  -  n  (2) 

1  -  n  (3) 

P  -  pa 

g  -  ga 

df 1  -  nfix  * (k  -1) 

df2  -  nfix  * (k* (1  -1)) 
irule  -  2 
c 

c  —  Find  peak  of  integrand  and  determine  transformation 

xO  -  sup  (float (20),  ier) 
if  (ier  .ne.  0)  xO  -  oldsup 


Number  of  fixed  factor  levels  (nfix) 
Number  of  batches  (k) 

Batch  size  (1) 

Number  of  intervals  for  k-function 
Values  of  nuisance  parameter 
Input  k-function 
Output  k-function 
Coverage  probability  of  xkl 
Probability  level  of  quantile 
Confidence  level  of  tolerance  limit 


A19 


if  (xO  .eq.  zero)  xO  -  dfl  / (df 1  +df2)  « 

oldsup  =  xO 

alpha  =  (dfl  /df2)  * (one  -xO)  /xO 
write  (*,*)  '  alpha  =  ',  alpha 
c 

do  10  i=l,  nstp 

eta  =  alpha*x(i)  tone 
c 

c  —  First  integral 

idrv  =  0 

call  dqdag  (xint,  tl,  th,  daerr,  drerr,  irule,  xiO,  errest) 
c 

c  —  Second  integral  (derivative) 

idrv  =  1 

call  dqdag  (xint,  tl,  th,  daerr,  drerr,  irule,  xil,  errest) 
cvrate  (i)  =  xiO 

xkstep  (i)  =  (g  -cvrate  (i) )  /xil 
xkl  (i)  =  xkO  (i)  +xkstep  (i) 

if  (diag  .eq.  one)  write  (*,*)  'nextk  :  k,  cvrate,  step  ' , 

$  i,  x  (i) ,  xkO  (i) ,  cvrate  (i) ,  xkstep  (i) 

10  continue 

c 

return 

end 

real  function  sup  (r,  ier) 

c  - 

c 

c  Mark  Vangel,  Dec  1988 

c 

c  Find  the  maximum  of  Trickett-Welch  integrand  for 

c  variance  ratio  equal  to  r.  The  spline  for  XK  must  be 

c  initialized  before  this  routine  may  be  used.  Also,  the 

c  stuff  in  /ca/  must  be  available, 

c 

common  /ca/  nfix,  k,  1,  eta,  xkp,  xkg,  idrv,  con 
external  fncs 

data  one  /!./,  eps  /l.e-5/,  xacc  /.001/ 
eta  =  l*r  +one 
c 

c  —  find  minimum  by  Brent's  method 

xguess  =  one  / 2 
bound  =  xguess  -eps 
xstep  =  one  / 4 
maxfn  =  100 

call  uvmif  (fncs,  xguess,  step,  bound,  xacc,  maxfn,  peak) 
sup  =  peak 
c 

return 

end 

subroutine  tlmcdf  (xmu,  s2b,  s2w,  n,  p,  t,  cum,  flenme) 

c  - 

c 

c  Mark  Vangel,  Oct.  1988 

c 

c  Cumulative  distribution  function  of  the  lower  confidence 

c  bound  on  s  quantile  from  a  random  effects  model.  This  routine 

c  can  be  used  with  a  nested  model.  The  results  in  the  paper 

c  correspond  to  nfix=l. 

c 

c  xmu  —  population  mean 

c  s2b  —  population  variance  between  groups 

c  s2w  —  population  variance  within  groups 

c  nfix  —  number  of  groups 

c  i  —  number  of  batches 
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C  j 

c  p 

c  t 

c  cum 

c  flenme 

c 
c 
c 
c 
c 
c 
c 

logical  known 
character*20  flenme 
dimension  n  (3) 
c 

c  —  Spline  for  tolerance  limit  factor 

common  /cb/  xs  (500),  xks  (500),  break  (500),  c  (4,  500),  nx 
c 

c  --  Flag  set  if  rho  taken  to  be  known 

common  /kw/  known,  dummy,  xknown 
data  lfn  /10/ 
c 

c  —  Initialize  spline  for  tolerance  limit  factor 

if  (.not.  known)  then 

open  (unit^lfn,  file=flenme,  iostat=lstat) 
rewind  (lfn) 

call  inispl  (lfn) 
end  if 

nfix  =  n  (1) 
i  =  n  (2) 
j  =  n  (3) 
c 

c  —  Set  up  parameters 

aerr  *  l.e-5 

xkp  =  anorin  (p) 

eta  »  j  *s2b  /s2w  +1. 

idf 1  =  nfix  *  (i-1) 

idf2  =  nfix  *(i  *(j  -1)) 
cl  =  (idfl  +idf2)  /nfix 
c2  =  cl  /eta 
cl  =  cl  *i  / (i  -1 . ) 

xncp  *  (xmu  -t)  /sqrt  ((j*s2b  +s2w)  /(i*j)) 
c 

c  --  Evaluate  cdf  of  lower  tolerance  limit 

call  evcdf  (cum,  idfl,  idf2,  cl,  c2,  xncp,  eta,  aerr) 
c 

return 

end 

real  function  welch  (ymsr,  xkg,  xkp,  m) 

c  - 

c 

c  Mark  Vangel,  July  1986 

c 

c  First  order  Welch-type  expansion  for  the  tolerance 

c  limit  factor, 

c 

real  i,  1 
dimension  m  (3) 
c 

i  -  m (1)  * (m{2) -1)  +1 
1  -  real  <m(l)  *m<2)  *(m(3)-l))  / 

$  real  (m(l)  *(m(2)-l)  +1)  +1 
c 


batch  size 

probability  associated  with  quantile 
value  at  which  cdf  is  to  be  evaluated 
probability  tol.  limit,  is  less  than  pth  quantile 
name  of  ASCII  file  containing  tolerance  limit  factor 
row  1  :  number  of  steps 

row  2..n  :  msr  / (msr  +1),  k-factor 
row  2  :  msr  ■=  0 
row  n  :  msr  —  infinity 
(output  from  porgram  TRICK) 
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<n  </></>  <a 


ymsr 


xmsr  = 
n  =  i*l 

tl  *  sqrt  (1/ (1+ (1-1) /xmsr) ) 

t2  =  sqrt  (1/ (xmsr*xmsr  + (1-1) *xmsr) ) 

rti  =  sqrt  (i) 

rtn  =  sqrt  (float (n) ) 

xll  -  1/  (1*1) 
xl2  =  (  (1-1) /l)  **2 
c 

xk  «=  xkp  +tl/rtn  *  (xkg  +1 .  /  (4*  < i— 1  >  )  *( 

xkg  *(xkg*xkg  +1)  +xkp*xkp*xkg  *n  *tl*tl  *xll 

+xkp  *rtn  *tl*tl*tl  *xll  +xkp*xkg*xkg  *rtn  *tl  /l) 

+1/ (4*i* (1-1 .dO) )  *( 

+xkp*xkp*  xkg  *n  *t2*t2  *xl2  +xkp  *rtn  *t2*t2*t2  *xl2) ) 
c 

welch  =  xk 

return 

end 

double  precision  function  xint  (x) 

c  - 

c 

c  Mark  Vangel,  Oct.  1988 

c 

c  Function  to  calculate  two  integrands  needed  for  the 

c  Trickett-Welch  procedure.  One  integrand  is  a  noncentral 

c  t  cumulative  'weighted'  by  a  beta  density;  the  other 

c  integrand  is  the  derivative  of  this  first  integrand  with 

c  respect  to  the  k-function  (which  is  part  of  the  argument 

c  of  the  noncentral  t  cumulative) . 

c 

double  precision  x 

common  /ca/  nfix,  k,  1,  eta,  xkp,  xkg,  idrv,  con 
data  one  / 1 .  / ,  half  /0.5/,  zero  /0.0/,  tiny  /l.e-6/ 
c 

c  —  Between,  within,  total  degrees  of  freedom 

df 1  =  nfix  * (k  -1) 

df 2  =  nfix  * (k  * (1  -1) ) 

df  -  dfl  +df 2 
c 

c  --  Calculate  mean  square  ratio.  Use  asymptote  when  mean 

c  square  ratio  is  infinite 

xl  =  x 

if  (xl  .le.  zero)  xl  =  tiny 
r  =  eta  *df2  *xl  /  (dfl  * (one  -xl) ) 

c 

c  —  Cubic  spline  interpolation  of  k-function 

call  kspln  (r,  xk) 
c 

c  --  Noncentrality  parameter  and  argument  for  noncentral  t 

xncp  =  xkp  *sqrt  (k*(one  +  (1  -one)  /eta)) 
arg  =  sqrt  (df  /nfix  *(k/(k  -one)  *xl  +(one  -xl)  /eta)) 

c 

c  —  This  subroutine  can  calculate  higher  derivatives  than 

c  the  first  if  desired, 

argl  =  one 

fact  “  one 

do  10  i»l,  idrv 

argl  *  argl  *arg 
fact  =  fact  *i 
argl  =  argl  /fact 
10  continue 

c 

c  --  Noncentral  t  cunulative  or  it's  derivatives 
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call  nctdrv  (idrv,  int (df ) ,  xk  *arg,  xncp,  prob) 


—  Beta  density 

if  (xl  *(one  -xl)  .ne.  zero)  then 
beta  =  (half  *dfl  -one)  *log  (xl) 

+(half  *df2  -one)  *log  (one  -xl) 
else  if  (xl  .eq.  zero)  then 

beta  =  (half  *df2  -one)  *log  (one  -xl) 
else 

beta  -  (half  *dfl  -one)  *log  (xl) 
end  if 

—  finally,  return  the  integrand, 
xint  =  argl  *prob  *exp  (beta  +con) 


return 
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